﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "ginv")]
    public static unsafe int ginv(IntPtr a_ptr, int m, int n, IntPtr aa_ptr, double eps, IntPtr u_ptr, IntPtr v_ptr, int ka, int inter = 100)
    {
        double* a = (double*)a_ptr.ToPointer();
        double* aa = (double*)aa_ptr.ToPointer();
        double* u = (double*)u_ptr.ToPointer();
        double* v = (double*)v_ptr.ToPointer();

        return ginv(a, m, n, aa, eps, u, v, ka, inter);
    }

    /// <summary>
    /// 求矩阵广义逆的奇异值分解法
    /// </summary>
    /// <param name="a">a[m]n]存放m×n的实矩阵A。返回时其对角线给出奇异值（以非递增次序排列），其余元素均为0。</param>
    /// <param name="m"></param>
    /// <param name="n"></param>
    /// <param name="aa">b[n][m]返回A的广义逆。</param>
    /// <param name="eps">给定的精度要求。</param>
    /// <param name="u">u[m][m]返回左奇异向量U。</param>
    /// <param name="v">v[n][n]返回右奇异向量V。</param>
    /// <param name="ka">其值为max(m，n)＋1。</param>
    /// <param name="inter">最多迭代次数，100</param>
    /// <returns>若小于0，则表示失败；若大于0，则表示正常。</returns>
    public static unsafe int ginv(double* a, int m, int n, double* aa, double eps, double* u, double* v, int ka, int inter = 100)
    {
        int i, j, k, l, t, p, q, f;

        i = muav(a, m, n, u, v, eps, ka, inter);
        if (i < 0) return (-1);
        j = n;
        if (m < n)
        {
            j = m;
        }
        j = j - 1;
        k = 0;
        //while ((k <= j) && (a[k * n + k] != 0.0))
        while ((k <= j) && (NE0(a[k * n + k])))
        {
            k = k + 1;
        }
        k = k - 1;
        for (i = 0; i <= n - 1; i++)
        {
            for (j = 0; j <= m - 1; j++)
            {
                t = i * m + j;
                aa[t] = 0.0;
                for (l = 0; l <= k; l++)
                {
                    f = l * n + i;
                    p = j * m + l;
                    q = l * n + l;
                    aa[t] = aa[t] + v[f] * u[p] / a[q];
                }
            }
        }
        return (1);
    }

#if __DRIVE_CODE__
      int main()
      { 
          int m,n,ka,i,j;
          double a[5][4]={ {1.0,2.0,3.0,4.0},
                      {6.0,7.0,8.0,9.0},{1.0,2.0,13.0,0.0},
                      {16.0,17.0,8.0,9.0},{2.0,4.0,3.0,4.0}};
          double aa[4][5],c[5][4],u[5][5],v[4][4];
          double eps;
          m=5; n=4; ka=6; eps=0.000001;
          cout <<"MAT A IS:\n";
          for (i=0; i<=4; i++)
          { 
              for (j=0; j<=3; j++)  cout <<a[i][j] <<"    ";
              cout <<endl;
          }
          i=ginv(&a[0][0],m,n,&aa[0][0],eps,&u[0][0],&v[0][0],ka);
          if (i<0)  return 0;
          cout <<"MAT A+ IS:\n";
          for (i=0; i<=3; i++)
          { 
              for (j=0; j<=4; j++)  cout <<aa[i][j] <<"    ";
              cout <<endl;
          }
          i=ginv(&aa[0][0],n,m,&c[0][0],eps,&v[0][0],&u[0][0],ka);
          if (i<0)  return 0;
          cout <<"MAT A++ IS:\n";
          for (i=0; i<=4; i++)
          { 
              for (j=0; j<=3; j++)  cout <<c[i][j] <<"    ";
              cout <<endl;
          }
          return 0;
      }
#endif
}
